function [res, param] = features(input, type, varargin)
% FEATURES generates various features from input data.
% Usage :
%   RES = FEATURES(INPUTS, TYPE, ...)
%   RES = FEATURES(INPUTS, TYPE, Param1, Val1, Param2, Val2, ...)
%
%   INPUTS is groups of input data in the form of cell.
%   Each input can be a matrix of image, a path to an image or a directory, 
%   or a search pattern. TYPE indicates types of features to extract. 
%   TYPE is case insensitive.
%   The following values of TYPE is supported:
%     'hist', 'point-hist', 'markov', 'diff-markov', 'abs-diff-markov',
%     'point-diff-markov', 'moment',
%
%   The following TYPE depends on third party tools:
%     'lbp*' depends on LBP Matlab code, which is available at
%   [www.cse.oulu.fi/CMV/Downloads/LBPMatlab].
%     'ccpev-274', 'ccpev-548' depends on CC-PEV feature extractor, which
%     is available at [dde.binghamton.edu/download/feature_extractors].
%
%   Other supported local texture features are:
%     LBP, LTP, LPQ, WLD, Gabor, IFLT
%
% Param and Values:
%   Arguments that can be passed to all TYPE are listed below:
%
%   Param Names | Value |      Description
% ------------------------------------------------------------------------
%  'use_jpegtbx'|   1   | Use jpegtbx to extract DCT from JPEG file
%                         if available. (by default)
%               |   0   | Do not use jpegtbx, always use imread/imwrite
%  'filetype'   |   *   | File extensions to scan in input directory.
%                         ('*.jpg' by default)
%  'qf'         | 1~100 | Quality factor for DCT. (75 by default)
%  'output'     |   *   | Output filename or type to store result.
%                         ('mat' for automatically generated names.)
%  'verbose'    |   2   | Print debug messages, processed files, and
%                         errors.
%               |   1   | Print processed files, messages and errors.
%                         (by default)
%               |   0   | Do not print any additional messages, except for
%                         errors.
%  'threshold'  |  num  | Threshold value or percentage. See also TYPE.
%  'points'     |  mat  | Points in DCT blocks to extract features. ([] by
%                         default)
%  'space'      | 'dct' | Input data in JPEG 8x8 DCT blocks.
%               | 'grey'| greyscale spatial image. (By default)
%               |'color'| Colorful/Multichannel spatial image.
%               | 'path'| Path to image files.
%               |'error'| Error prediction image.
%  'diff_dir'   |'x'|'y'| Difference in horizontal/vertical direction.
%  'diff_dis'   |  >=1  | Difference distance.(1 by default)
%  'diff_n'     |  >=0  | Difference times.(0 by default)
%  'pe_func'    | @???? | Callback for error image calculation. Default 
%                         callback is one proposed in papers by Shi et al.
%  'pe_n'       |  >=0  | Error prediction times.(0 by default)
%  'use_abs'    | 1 / 0 | Use absolute value instead of original data.
%                         (0 by default)
%   Preprocess is in this order:
% ... -> absolute value -> difference -> prediction error -> ...
% Type-specific information:
%   'hist' :
%     generates histogram of all coeffients. Supported param
%     names are listed below:
%
%     Param Names |       Description               |  Default Value
%     -----------------------------------------------------------------
%     'threshold' | Threshold value or percentage.  |  10
%     'space'     | 'dct' or 'grey'                 |  'dct'
%
%     If threshold is a positive number, T is equal to threshold.
%     Or if it is a float value between 0..1, T is determined that
%     at least threshold% fall in range.(TODO)
%     Output feature is a (., 2*T+1) matrix, standing for probability
%     of values in range [-T:+T].
%
%   'point-hist' :
%     generates histogram of each coeffient in DCT. Supported param
%     names are listed below:
%
%     Param Names |       Description               |  Default Value
%     -----------------------------------------------------------------
%     'threshold' | Threshold value or percentage.  |  10
%     'points'    | Points in DCT block.            |  []
%
%     If points is an 2*N matrix. The first line stands for ROW, the
%     second line stands for COLUMN, and the value in it is in range
%     1..8. Otherwise, if points is an 1*N matrix. The points is counted
%     from 0..63, 0 stands for DC, 1..63 stands for AC, in zig-zag order.
%     Output feature is (N, .) matrix, one line of output corresponds
%     to one point in points.
%
%   'diff-markov' :
%    Directional Difference Markov Feature.
%                    Output feature is in range [-T:T, -T:T].
%
%   'abs-diff-markov':
%     generates directional difference Markov feature. It is only
%     available if space is DCT. It uses absolute value of DCT
%     coeffients.
%
%     Param Names |       Description               |  Default Value
%     -----------------------------------------------------------------
%     'threshold' | Threshold value or percentage.  |  4
%
%     Output feature is in range [-T:T, -T:T].
%
%   'point-diff-markov' :
%     Not implemented.
%
%   'moment' :
%     calculates moments of input.
%
%     Param Names |       Description               |  Default Value
%     -----------------------------------------------------------------
%     'space'     | 'dct' or 'color'                |  'dct'
%
%   'dlbp' :
%     Param Names |       Description               |  Default Value
%     -----------------------------------------------------------------
%     'ratio'     | Ratio of Dominant patterns      | 0.8
%     'k_ratio'   | Average value of K              | []
%     'radius'
%     'neighbors'
%     'ntrain'
%
% Return value:
%   RES is M*N matrix, where N is the number of features and M is the
%   number of samples.
%   PARAM is a struct. It contains input parameters.
%
% Change logs:
% r26, 04/21/12, initial version, support histogram and markov.
% r37, 05/14/12, use varargin and struct to parse/save inputs.
% r44, 06/30/12, add point histogram, point difference markov.
% r72, 07/24/12, add libsvm output, add lbp, cc-pev, moment features.
% r95, 08/24/12, add gabor, lpq, wld, ltp, iflt features.
% r97, 09/15/12, add predict error image.
% r124, 03/25/13
%  * add diff_dir/diff_n/diff_dis/pe_n/pe_func.
%  * return initial settings now.
% current, 03/26/13, code & doc cleanup, adding dlbp
%
% AUTOCOMPLETION:
%   Run "edit(fullfile(matlabroot,'toolbox/local/TC.xml'))" and append following lines:
%
%     <binding name="features">
%         <arg argn="2" ctype="VAR" value="hist blockhist markov markov2-noncausal diff-markov abs-diff-markov moment ccpev-274 ccpev-548 lbp gabor lpq ltp wld iflt"/>
%     </binding>

%% parse input arguments
[param, type] = default_args(type);
param = parse_args(param, varargin);
param.type = type;

if ischar(input) % path
    if exist(input, 'file') && (~exist(input, 'dir')) % file
        if param.verbose > 1
            fprintf(1, 'Treat as file: %s.',input);
        end
        inlist = dir(input);
        indir = fileparts(input);
    elseif exist(input, 'dir') % directory
        if param.verbose > 1
            fprintf(1, 'Treat as directory: %s.', input);
        end
        inlist = dir([input,'/', param.input_filetype]);
        indir = input;
    else % pattern?
        if param.verbose > 1
            fprintf(1, 'Treat as pattern: %s.', input);
        end
        inlist = dir(input);
        indir = fileparts(input);
    end
    
    M = size(inlist, 1);
    if M == 0
        warning('features:NoMatchedFile','No such files or directory: %s.', input);
    end
    if param.verbose > 1
        fprintf(1, 'Found %d matches for %s.',M, input);
    end
    
    switch lower(type)
        case 'hist'
            if isfield(param, 'threshold') && (param.threshold < 1)
                param = pre_scan(param, type, inlist, indir);
            end
        case 'point-hist'
            if isfield(param, 'threshold') && (param.threshold < 1)
                param = pre_scan(param, type, inlist, indir);
            end
        case 'markov'
            if isfield(param, 'threshold') && (param.threshold < 1)
                param = pre_scan(param, type, inlist, indir);
            end
        case 'dlbp'
            if ~isfield(param, 'k_ratio') 
                param = pre_scan(param, type, inlist, indir);
            end
        otherwise
            %pass
    end
    
    param.feature_size = get_feature_size(type, param);
    if param.feature_size > 0
        % known feature size
        res = zeros(M, param.feature_size);
    end
    param.global = prepare_global(type, param);
    if param.verbose > 0
        disp(param);
    end
    count = 0;
    
    for i = 1:M
        if (inlist(i).isdir ~= 0) || (inlist(i).bytes <= 0)
            continue;
        end
        name = inlist(i).name;
        fullname = [indir,'/',name];
        [v, succ] = deal_one(fullname, type, param);
        if succ
            if param.feature_size <= 0
                % decide feature size in first run
                param.feature_size = numel(v);
                res = zeros(M, param.feature_size);
                disp(param.feature_size);
            end
            count = count + 1;
            res(count, :) = v(:)';
            if param.verbose > 0
                fprintf(1, 'Processed %s.\n', name);
            end
        else
            if param.verbose > 0
                fprintf(1, 'Failed %s.\n', name);
            end
        end
    end
    res = res(1:count,:);
else
    % matrix?
    param.feature_size = get_feature_size(type, param);
    if param.feature_size > 0
        % known feature size
        res = zeros(1, param.feature_size);
    end
    if param.verbose > 0
        disp(param);
    end
    res = feature_extractor(input, type, param);
end
fprintf(1, 'All input processed.\nSave results:\n');

% output file
switch lower(param.output)
    case ''
        % pass
    case 'mat'
        param.output = sprintf('%s/%s_%d.mat', input, type, param.threshold);
        if isfield(param,'output_append') && (param.output_append ~= 0)
            save(param.output, '-append', 'res');
        else
            save(param.output, 'res');
        end
        fprintf(1, 'saved to %s\n', param.output);
    case 'libsvm'
        param.output = sprintf('%s/%s_%d.dat', input, type, param.threshold);
        if isfield(param,'output_append') && (param.output_append ~= 0)
            fout = fopen(param.output, 'a');
        else
            fout = fopen(param.output, 'w');
        end
        if fout == -1
            error('features:output:libsvm:CannotWriteFile', 'Cannot write output file %s.', param.output);
        end
        [M,N] = size(res);
        for i = 1:M
            fprintf(fout, '%d', param.outputclass);
            for j = 1:N
                fprintf(fout,' %d:%0.6f',j,res(i,j));
            end
            fprintf(fout,'\n');
        end
        fclose(fout);
    otherwise
        save(param.output, 'res');
        fprintf(1, 'saved to %s\n', param.output);
end
end %features

% param is arguments passed in
% type is feature type
% inlist is file list
function param = pre_scan(param, type, inlist, indir)
fprintf(1, 'Prepare pass: scan files and determine stats.\n');
count = 0;
% Defining callback function to prepare values:
%
% _tmp = XXX_prepare_start(param)
%   prepare data for processing
%
% _tmp = XX_prepare_process(img, param, _tmp, count)
%   process images one by one, all data need to pass must be returned
%
% param = XX_prepare_finish(param, _tmp, count)
%   store results into param, it is then returned
%
startfunc = 0;
processfunc = 0;
finishfunc = 0;
switch lower(type)
    case 'dlbp'
        startfunc = @DLBP_prepare_start;
        processfunc = @DLBP_prepare_process;
        finishfunc = @DLBP_prepare_finish;
    case 'hist'
        error('pre_scan not implemented');
        %% TODO: determine threshold
    case 'point-hist'
        error('pre_scan not implemented');
        %% TODO: determine threshold
    case 'markov'
        error('pre_scan not implemented');
        %% TODO: determine threshold
    otherwise
        warning('pre_scan should not be called for type %s', type);
        return;
end
M = length(inlist);
perm = randperm(M);

if isa(startfunc, 'function_handle')
    data = startfunc(param);
else
    data = struct();
end

for i = 1:M
    if (inlist(perm(i)).isdir ~= 0) || (inlist(perm(i)).bytes <= 0)
        continue;
    end
    name = inlist(perm(i)).name;
    disp(name);
    fullname = [indir,'/',name];
    try
        img = imread(fullname);
        if isa(processfunc, 'function_handle')
            data = processfunc(img, param, data, count);
        end
        count = count + 1;
    catch ME
        ME
    end
    % scan total histogram of all files
    % returned: 1)count 2)param.threshold
    fclose('all');
end

if isa(finishfunc, 'function_handle')
    param = finishfunc(param, data, count);
end

if param.verbose > 0
    fprintf(1, '%d files scaned.\n', count);
end
end